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This invention generally relates to geologic fault surfaces. More specifically, 
this invention is a method to interpret and construct a geologic fault surface to fit 
three-dimensional seismic discontinuity data. 



In hydrocarbon system evaluations, mapping faults and fault networks is 
essential to determine the migration pathways from the source to the reservoir. Faults 
can also help trap hydrocarbons or fragment a reservoir and therefore cause 
complications during field production. 

Fault interpretation and fault network interpretation in three-dimensional 
seismic data can be facilitated and accelerated by the use of seismic amplitude 
discontinuity data. For example, U.S. Patent Application No. 09/827,574 (Cheng et 
al.) discloses a method of identifying structural and stratigraphic discontinuities in a 
three-dimensional seismic data volume through the use of seismic amplitude 
discontinuity data. In Cheng et al, the continuity of seismic reflectors in a volume of 
seismic amplitude data is measured by computing the correlation coefficient between 
adjacent seismic traces over a movable vertical window. A low coefficient of 
correlation indicates that the reflector is discontinuous. Repeating this practice over an 
entire volume of seismic data creates a discontinuity cube characterizing the 
continuity of the reflectors in the seismic volume. Since faults are detected by looking 
at vertical offsets of seismic reflectors, the discontinuity cube is a preferred way to 
image faults and fault networks in a volume of seismic data. 

Other methods for detecting faults include converting seismic data volumes 
into cubes of coherency measurements herein referred to as "seismic coherency 
cubes." This method is disclosed by U.S. Patent Nos. 5,563,949 and 5,838,564 
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(Bahorich and Farmer), which are commonly known as the "coherency cube" patents. 
For purposes of this application, seismic coherency cube data and seismic 
discontinuity cube data can be used interchangeably as fault-indicating parameters in 
the inventive method claimed in this application, both as a substitute for each other or 
5 in combination. 

Current technology in fault interpretation focuses on automatic fault detection 
and extraction using amplitude and coherency data. One characteristic of automatic 
fault detection is that no preexisting fault interpretations are required for automatic 
fault detection. However, a key issue in automatic fault detection is the quality of the 

10 seismic amplitude and coherency data used to detect the faults. Therefore, automatic 
fault detection methods require preprocessing of seismic data to enhance the quality of 
the fault signature in amplitude and coherency data and to facilitate generation of 
specific criteria to facilitate differentiation of true faults from false fault signatures 
during the extraction process. For example, U.S. Patent No. 5,987,388 (Crawford et 

15 al.) discloses an automatic fault detection method. Another approach based on 
mathematically inserting test planes into a volume of seismic data to approximate dip 
and azimuth of potential fault surfaces is disclosed in U.S. Patent No. 6,018,498 (Neff 
et al.). 

Automatic fault detection may work well with extremely good quality data. 
20 However, many seismic data volumes do not have the quality required for automatic 
methods. Therefore, auto-assisted methods where the seismic interpreter guides the 
computer by inserting partial interpretations generally are more reliable, particularly 
with data of lesser quality. 

US Patent No. 5,537,320 (Simpson et al.) discloses one example of an auto- 
25 assisted method. This method starts with a manually interpreted fault stick line (a 
piecewise linear line), which is defined by at least two fault nodes, in a vertical slice 
of a seismic amplitude volume. A "fault stick** represents the intersection of a fault 
surface and any planar slice of the data volume. The nodes are points in the slice 
lying along the fault, and are typically identified by a seismic interpreter. The method 
30 of Simpson et al. requires several processes that are used to refine and extend the 
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initial fault seed nodes or initial points designated by the user to represent the fault. 
First, a "snap" process moves the fault seed nodes so as to be located at voxel points 
at which minimum correlation occurs between seismic amplitudes in either side of the 
fault nodes. Voxel points are points in space (or in a grid of a three-dimensional 
5 volume) with a location (such as (x,y,z) coordinate) and value (typically, grayscale 
from 0-255) representing seismic amplitude or its discontinuity. The next step is to 
extend the two end-point fault seed nodes. An end-point fault seed node is a fault 
point located at the two ends of a fault stick (or fault polyline). The two end-point 
fault seed nodes are extended in upward and downward directions respectively with a 
10 fixed length and the process makes a decision if the two end-nodes can be extended. 
The final step projects the snapped and extended fault nodes to a next vertical slice. 
These projected fault nodes serve as new fault seed nodes and the process is repeated. 
In this three-step process, a quality control threshold value is used to stop extensions 
in vertical directions and into the next slices. 

15 The Simpson et al. processes do not use the three-dimensional information 

inherent in the fault surfaces in a seismic amplitude cube. The "snap'* process is a 
two-dimensional operator, meaning that the decision on moving a fault node in one 
vertical slice is based only on the information obtainable from that vertical slice. 
Even in a given vertical slice, the movement of one node is not constrained in any way 

20 by the location of the other nodes in the same vertical slice. Furthermore, the 
technique is not able to handle the fault nodes selected in horizontal slices and vertical 
slices together. Segmentation of the shapes and boundaries of a three-dimensional 
object based on a two-dimensional operator without sufficient smoothness constraints 
is known in the art to be very susceptible to noise in voxel values. 

25 In one variation, disclosed in Simpson et al., fault seed nodes in two or more 

vertical slices are jointly interpolated to generate fault seed nodes for intervening 
vertical slices. These newly created nodes in each vertical slice are refined by using 
the "snap" process in each vertical slice. 

The result of these prior art hydrocarbon system oriented techniques is the 
30 characterization of the object of interest, generally a fault surface, in the three- 
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dimensional data set. The problem of finding and parameterizing shapes and 
boundaries of an object in two- and three-dimensional images has also been 
extensively studied in the image analysis and computer vision literature. Promising 
models that have robustness to noise and the flexibility to represent a broad class of 
5 shapes include deformable surface models and their two-dimensional analog, active 
contours. See, for example, M. Kass, et al. "Snakes: active contour model," Int. J. 
Comput. Vision 1, 321-331 (1988). Kass discloses two-dimensional models, but the 
principles for three-dimensional models are the same. A deformable surface behaves 
like an elastic sheet. Initially placed close to a boundary of interest, a deformable 

10 surface deforms towards the desired boundary under the influence of external forces 
(attraction toward salient voxel features or anomalies in the data such as high 
discontinuity) and internal elastic forces (surface smoothness constraints). Variations 
of deformable models have been successfully utilized in reconstructing boundaries of 
brain, heart tissue, and blood vessels from medical images. Examples of references 

15 disclosing three-dimensional models for medical applications include L. Cohen and I. 
Cohen, "Finite-element methods for active contour models and balloons for 2-D and 
3-D images," IEEE Trans. PAMI 15, No. 1 1 (1993); A. Bosnjak, et al., "Segmentation 
and VRML visualization of left ventricle in echocardiograph images using 3D 
deformable models and superquadrics," Proceedings of the 22 nd Annual International 

20 Conference of the IEEE Engineering in Medicine and Biology Society 3, 1724-1727 
(July, 2000); and M. Hernandez-Hoyos et al., "A deformable vessel model with single 
point initialization for segmentation, quantification and visualization of blood vessels 
in 3D MRA," Proceedings of the Third International Conference of Medical Imaging 
Computing and Computer-Assisted Intervention, 735-745 (Oct. 2000). 

25 Patents have been issued for using two-dimensional deformable models for 

automatically determining the boundaries of objects in three-dimensional topographic 
images (U.S. Patent No. 6,249,594), and using wavelet coefficients and a two- 
dimensional model for the detection of nodules in biological tissues (U.S. Patent No. 
6,078,680). Deformable models have been applied to geoscience by Apprato et al, 

30 who used a traditional two-dimensional deformable model on a sea floor image to 
detect a fault line. D. Apprato, C. Gout, and S. Vieira-Teste, "Segmentation of 
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complex geophysical 3D structures," Proceedings of IEEE 2000 National Geoscience 
and Remote Sensing Symposium, p. 651-653, July, 2000. Apprato et al uses a 
standard two-dimensional deformable model for fault segmentation, but on a two- 
dimensional sea floor image. The model is not suitable for a three-dimensional fault 
5 surface segmentation because the model does not utilize seismic volume data that is 
required for three-dimensional fault surface construction. 

Due to the above-mentioned difficulties with traditional automatic 
interpretations, faults are typically interpreted manually using both amplitude and 
discontinuity data. The interpreter scrolls through the vertical section of the seismic 

10 amplitude cube and digitizes the fault sticks. The fault interpretation is simultaneously 
co-referenced on time slices of the discontinuity cube for verification that the 
interpretation satisfies the fault trace on the discontinuity time slice. Depending on the 
level of accuracy required and the complexity of the fault network, the operator or 
interpreter may increase the number of fault sticks necessary to describe the fault. 

15 Finally, the interpreted fault is gridded or triangulated to create a fault surface using 
commercially available software (i.e., Gocad™ or Voxelgeo™). The problem with 
existing methods is that traditional manual interpretation methods require time 
consuming dense fault sticks for an accurate fault surface and automatic interpretation 
methods lack accuracy in three-dimensions. Accordingly, there is a need for a rapid, 

20 accurate, fault interpretation method. The present invention satisfies that need. 



SUMMARY 



A method to create fault surfaces from a three-dimensional seismic data 
volume is disclosed. In one embodiment, the method comprises: (a) generating at 

25 least two fault sticks containing interpreter-provided fault nodes, each fault stick 
coming from a different slice of the data volume; (b) constructing an initial three- 
dimensional fault surface that includes the fault sticks; and (c) reconstructing the fault 
surface to fit the seismic amplitude discontinuity or coherency information in the data 
by iterating a deformable surface model to obtain a realistically smooth surface that 

30 passes through voxels having high values of seismic amplitude discontinuity or low 
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values of amplitude coherence. In some embodiments, the iteration of the deformable 
surface model is accomplished by local minimization of an energy function of the 
surface, where the surface energy is a function of the curvature of the fault surface and 
either the degree of seismic amplitude discontinuity or the degree of seismic 
5 amplitude coherency on the fault surface. In some embodiments, the surface energy 
function is a weighted sum of an internal force being a function of surface curvature 
and representing a smoothness constraint, and an external force being a function of 
discontinuity or coherency and tending to cause the surface to pass through voxels 
having high values of discontinuity or low values of coherency. Specific energy 
10 functions in the form of weighted sums are disclosed as examples. A preferred 
method of constructing the initial fault surface from the user-provided fault sticks is 
disclosed. A method of speeding up the iteration of the deformable surface model is 
disclosed taking advantage of the fact that fault surfaces tend toward being vertical 
planar surfaces. 

15 BRIEF DESCRIPTION OF THE DRAWINGS 

Fig. 1 shows a flowchart of an embodiment of the present invention for fault 
surface construction; 

Figs. 2A and 2B show a three-dimensional fault surface represented by y = f(x, 

z); 

20 Fig. 3A shows a fault surface in the x-y plane; 

Fig. 3B shows a fault surface in the y-z plane; 

Fig. 4A shows surface curvature of a fault surface in the y-z plane; 

Fig. 4B shows surface curvature of a fault surface in the x-y plane; 

Fig. 5 shows a target fault surface and an initial fault surface; 

25 Fig. 6A shows a horizontal slice of the target fault surface and initial fault 

surface; 
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Fig. 6B shows a horizontal slice of the initial fault surface and the final fault 
surface overlapped over the target fault surface; 

Figs. 7A, 7B, 7C and 7D show a horizontal slice of the target and initial fault 
surfaces and the movement of the constructed fault surfaces at different stages of the 
5 evolution when 30% noise is added to the discontinuity voxel volume; 

Figs. 8A, 8B, 8C, and 8D show a horizontal slice of the target and initial fault 
and the movement of the constructed fault surfaces when 50% of the target fault is 
missing and 20% noise is added to the discontinuity voxel volume; 

Figs. 9 A, 9B, 9C, and 9D show a graphical flowchart of the steps of an 
10 embodiment of the invention using publicly available three-dimensional seismic data 
from the Gulf of Mexico. 

DETAILED DESCRIPTION 

In the following detailed description, the invention will be described in 
connection with its preferred embodiments. However, to the extent that the following 
15 description is specific to a particular embodiment or a particular use of the invention, 
this is intended to be illustrative only. Accordingly, the invention is not limited to the 
specific embodiments described below, but rather, the invention includes all 
alternatives, modifications, and equivalents falling within the true scope of the 
appended claims. 

20 One embodiment of the method disclosed below comprises four steps. First, at 

least two fault interpretations are generated. Preferably, the fault interpretations are 
generated manually using widely spaced or sparse fault sticks. Second, an initial fault 
surface is constructed. The fault surface may be generated automatically by bilinearly 
interpolating the fault sticks. Third, a fault surface is reconstructed to fit the seismic 

25 amplitude discontinuity data. Preferably, in this step, the initial fault surface is 
reconstructed to fit the discontinuity data through use of a modified deformable 
surface algorithm. Fourth, fault sticks or surface points are outputted (such as 
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displayed) at user or operator defined intervals. This method creates fault 
interpretations faster and with higher accuracy. 

The first step requires the operator to manually generate fault interpretations 
using fault sticks. The operator can generate the fault interpretation manually using 
5 techniques known to persons skilled in the art. In one embodiment of this invention, 
the fault surfaces are generated (the second step in the previous paragraph) by using 
operator provided sparse fault interpretations in vertical as well as horizontal slices of 
a seismic volume as input. The inputted interpretations are then used to define the 
areas of the fault surface and construct an initial fault surface. This initial surface is 
10 subsequently improved or reconstructed by minimizing a cost function that is 
composed of the smoothness of the surface and the degree of discontinuity in the 
surface. 

The fault surface is preferably reconstructed (the third step as discussed 
previously) utilizing a deformable surface model. The deformable surface model is 
15 specifically formulated for constructing fault surfaces, which are substantially vertical 
and smooth without meandering complex curvatures. 

The philosophy of the deformable model approach is to introduce an initial 
fault surface in the vicinity of a desired fault surface and let it evolve from an initial 
position. The evolution is under the action of both internal forces (smoothness 

20 constraints) and external forces (attraction toward salient voxel features such as high 
discontinuity). Therefore, the energy function of a surface is formulated to represent 
these two interacting forces. Then the evolution or the deformation of a surface takes 
place as the energy function is minimized until it reaches a local minimum energy 
state. Since the optimization guarantees only a local minimum, an initial fault surface 

25 is preferably located fairly near the desired fault boundaries for an effective fault 
surface reconstruction. 

The objective of using this deformable model is to construct a "smooth 
surface" that passes through high discontinuity voxels in a given region of 
discontinuity voxel volume. Compared to the Simpson et al. process, this method will 
30 be more robust to noise and produce more detailed fault surfaces. The simultaneous 
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optimization of the local curvatures and the degree of discontinuity in the fault surface 
in three-dimensional space provides for more detailed fault surfaces and additional 
robustness to noise. Compared to other deformable models, the method disclosed by 
this invention is faster with higher resolution by taking advantage of the generally 
5 smooth and vertical characteristics of fault surfaces. 

Fig. 1 is a flow chart illustration of this computer aided fault surface 
construction process. First, seismic discontinuity data volume 1 and seismic 
amplitude data volume 3 are inputted into the computer. The operator manually picks 
fault surface points 5 from the inputted data (such as, seismic amplitude data volume 
10 or discontinuity data) by drawing fault sticks from the data inputted in the first step. 
Next, an initial surface is generated through interpolation 7. Finally, a fault surface is 
reconstructed to fit the seismic data. The reconstruction is preferably generated 
through an iterative evolution of at least one fault surface through energy 
minimization 9. 

15 As discussed above seismic data is inputted into the computer. As shown in 

Fig. 1, the data consists of seismic amplitude data volume 3 and seismic discontinuity 
data volume 1. Furthermore, and as discussed above, coherency data volume can be 
used as a replacement for discontinuity data volume for purposes of this invention. 
For this invention there is no preprocessing required. However, in a preferred mode, 

20 seismic data can be preprocessed by a special noise reduction filter (such as a median 
filter) to improve the discontinuity cube data. Persons skilled in the art will recognize 
known processing techniques for improving the seismic data. 

As shown in Fig. 1, the operator manually picks fault sticks in at least two 
slices 5. The slices can be either vertical or horizontal slices and can be picked at 
25 random or the slice can be picked to represent features of interest. Preferably, the 
slices picked by the operator would be at least one horizontal slice (x-y slice) and at 
least two vertical slices (x-z and y-z slices). More preferably, the slices would 
represent at least one x-y slice, at least one x-z slice, and at least one y-z slice. 
However, the slices can be two slices of parallel planes. Increasing the number of 
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slices with manually picked fault surfaces would increase the number of fault features 
represented and would increase the accuracy of the method. 

A set of fault points (or fault nodes) that belong to an initial fault surface are 
provided by a manual fault surface point picking process in vertical as well as 
5 horizontal slices of a seismic volume. In step 7 of Fig. 1 , these initial fault surface 
points are interpolated to generate an initial fault surface 7. Signatures of fault surface 
points are best-detected in horizontal and vertical slices of a seismic amplitude and 
discontinuity data volume. Interpreters may use any graphical interface tool (such as 
Voxelgeo™) to pick fault surface points at vertical and/or horizontal slices and store 
10 these surface point coordinates into a computer's memory. The initial fault surface 
points from step 5 of Fig. 1 are picked from two-dimensional slices but the fault 
surface resulting from step 7 is generated through interpolation in three-dimensions. 

In most three-dimensional deformable models a fault surface is defined by a 
set of discrete points. These discrete points are later moved to the desired boundaries 

15 of an object by minimizing surface energy, which is a function of image pixel values 
and curvature at the surface. In the prior art processes, considerable computation time 
is required to resample the discrete points because the distribution of discrete points 
becomes irregular as they are moved in normal directions to surface curvatures. 
Furthermore, the computation of curvatures becomes very complex and time 

20 consuming. In some deformable models, finite element methods are used to reduce the 
complexity and the computing time of the curvatures. 

Fig. 2A illustrates a data volume containing a fault surface Q defined at 
regular intervals inside the domain R in the x-z plane. The three axes of the data 
volume are represented by the x-y-z arrows 25. The fault surface Q is bounded by 
25 points A, B, C and D. The dotted line connects the fault surface with domain R to 
illustrate the projection of the fault surface points A, B, C, and D onto the domain R 
bounded by points A', B', C and D\ For illustrative purposes, the value of x 
increases from A to B and x increases from D to C on the fault surface Q . 

In order to take advantage of the relatively simple shape of geologic faults 
30 (such as their substantially vertical, smooth, and non-closing surfaces), a fault surface 
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in this invention is assumed to be representable by a single valued function in the form 
of x =f(y, z) oxy =f (x, z) 9 where x andy are geographic coordinates in a horizontal 
plane and z is a vertical depth (or time) coordinate. A functional form of y =f(x, z) 



will be used if the angle 0 between the intersection of the fault surface Q and a plane 
5 parallel to the y-z plane is greater than 45 degrees. Otherwise, a functional form x =/ 
(y, z) will be used to represent a fault surface. Fig 2B illustrates a fault surface Q 
represented by y = / (x, z), where the angle 6 between the intersection of the fault 
surface Q and a plane parallel to the y-z plane 23 is assumed to be greater than 45 
degrees. In this invention, as shown in Fig. 2A, a fault surface function y =f (x, z) 
10 will be defined at regular intervals inside the domain R in the x-z plane, which is a 
projection of a fault surface Q onto an x-z plane. 

Fig. 3A and Fig. 3B illustrate resulting discrete voxel representations of a fault 
surface 21 function y =/ (x, z) at a horizontal slice (x-y plane) and a vertical slice (x-z 
plane, or y-z plane) respectively. In this example, a y-z plane was used to represent a 

15 vertical slice. It should be noted that, at both slices, y is a single valued function of x 
and z. With this surface representation, all the voxels in a fault surface will be allowed 
to move only along a y direction (x direction for a fault surface with x = / (y, z) 
representation). In other deformable models (such as, U.S. Patents Nos. 6,078,680 
and 6,249,594), sparsely located surface points are allowed to move in normal 

20 directions to the surface and the mathematical function describing a surface becomes 
more complex than the one in this invention. This simplified approach takes 
advantage of the relatively simple shape of geologic faults and eliminates the need for 
a complex surface representation, complex curvature computation, and re-sampling of 
fault surfaces at each iteration of the energy minimization process. 

25 Given a set of initial fault surface points, the functional form of a fault surface 

(either x = / (y, z) or y = / (x, z)) is determined by estimating the angle 6 between a 
best fitting plane through the fault surface points and y-z plane. For simplicity, but not 
meant to be limiting, the example illustrated in Fig. 2B, has an angle 9 greater than 
45 degrees and uses a fault surface function y =f(x t z) in describing the invention. 
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An initial fault surface can be generated (7 in Fig. 1) in three steps. In the first 
step, a linear interpolation is determined between initial user provided fault surface 
points in at least two slices (preferably, at least one horizontal and at least two vertical 
slices). Second, the extent of the fault surface boundaries are identified by using the 
5 interpolated fault points in the previous step. For the example (illustrated in Fig. 2A) 
of fault surface function y = / (x, z) 9 this process pertains to the definition of the 
domain R in (x, z). Finally, an initial fault surface is constructed by using bilinear 
interpolation of interpolated fault points in the previous step (step 5 of Fig. 1). 

As shown as step 9 in Fig. 1, the final step is a reconstruction of the fault 
10 surface to fit the seismic data. Preferably, the reconstruction is through an iterative 
evolution of the fault surface through energy minimization. The energy minimization 
is preferably accomplished through use of a deformable model. 

In the present example, a deformable surface model is used for the iterative 
evolution of a fault surface through energy minimization as shown as the final step 9 
1 5 in Fig. 1 . The use of a deformable surface model requires an initial fault surface to be 
located in the vicinity of an unknown fault surface. Then, the initial fault surface is 
deformed toward a final fault surface by minimizing the energy of the fault surface in 
an iterative way 9. In this example, the energy of a fault surface is defined as a 
weighted sum of the discontinuity voxel values, which represent the degree of 
20 discontinuity in seismic amplitudes at a fault surface, and the curvatures of the fault 
surface at voxel locations belonging to the fault. 

One embodiment of an energy function of a surface Q may be formulated as a 
weighted sum of internal force v s and external force p s defined at points s in regular 
intervals of (x, z) in Q : 

25 E= v s + 0-A) p s ) 0) 

seO 

where X is a weight between the two forces v s and p s . The internal force v s at a point 
5 in a surface Q represents smoothness constraint. The smoothness constraint 
emulates elastic property of a surface at a point s. Thus, limiting the extreme 
curvature of the surface. In this embodiment, v s is defined as the sum of the 
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smoothness measured in vertical and horizontal directions. Let s be the point at (jc, y, 
z). Then v s is defined as: 

v s = + v sh (2a) 
v J V = 1-cos 0 s v (2b) 
5 = 1-cos^ (2c) 

where, as depicted in Fig. 4A, 0 S V is an angle between two vectors connecting point 
41 with coordinates (x, y u , z- / v ) with point 43 with coordinates (x, y, z) and point 45 
with coordinates (x, j rf ,z+/ v ) inanx slice. In the present example, y u and y d are 
the y values of the single valued fault surface function y =f(x, z) at (x, z-l v ) and (x, 

10 z+/J respectively. As shown in Fig. 4B, the angle 0 jA is an angle between two 
vectors connecting points (x-/ A , y r , z) 47, (x, y, z) 43, and (x+/ A , y /9 z) 49 in a z 
slice. Again, y r and y t are the corresponding values from a single valued fault 
surface function j; = /(;c,z)at (x-l h ,z) and (x + l h9 z) respectively. The parameters 
l v and / A are chosen by the operator and represent respectively the vertical and 

15 horizontal lengths of a window, for which internal force and external force are 
computed for each fault surface point s. This internal force serves to impose 
smoothness on the fault surface. 

The external forces pushing the surface toward salient voxel features such as 
low coherency (inverse if discontinuity data is utilized) is defined as: 

20 Ps = Ps, v + Ps* (3) 

where, as shown in Fig. 4A, p s v is the sum of coherency (inverse of discontinuity) 
values along the lines connecting point 41 with coordinates (x, y u , z-/ v ) with point 
43 with coordinates (x, y, z) and point 45 with coordinates (x, y d9 z+/ v ). 
Furthermore, as shown in Fig. 4B, p s h is the sum of coherency voxel values along the 
25 lines connecting points (x-/ A , y r , z) 47, (x, y, z) 43, and (x+ l h , y, , z) 49 in a z slice. 
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At the beginning of a fault surface construction, the energy function E 0 of an 
initial fault surface Q 0 is computed according to Eq. (1). Then, every fault surface 
point s in Q 0 is allowed to move up to ± d voxel units in a direction parallel to the y 
coordinate. Among the (2 d+l) possible voxel locations, the best move for the surface 
5 point s that minimizes the energy function is selected as d* . Then, the next fault 
surface is formed by moving each surface point by preferably an a fraction of d s , 
where a is a minimization step size. A value of a less than 1 promotes a stable 
minimization process. 

This iterative minimization repeats until the fault surface energy does not 
10 change any more or the change is less than a predetermined value. In another 
embodiment, one can smooth the voxel movement by applying weighted averages of 
the optimum voxel movements. Instead of moving the point s by a d\ , the point s is 
moved by a d'* , where d? is a weighted average of best moves of the surface points 
in the vicinity of the point s and weights are inversely proportional to the distances 
15 between the point s and the neighboring points. The smoothing window, which 
defines the neighboring points in domain R, as shown in Fig. 2 is denned as a 
rectangle that is centered at s with a width 2 r h +1 and a length 2 r v +1 . Even though the 
size of this smoothing window could differ from the window size for computing 
internal and external forces, the same size: r h =l„ and r v =/„ is generally acceptable. 

20 The weighting factor X in Eq. (1) determines the elastic property of the 

surface. The higher the value the more rigid the surface becomes. An extreme value 
of X =1 will force the surface to become close to a flat plane. A low value of X will 
let the surface follow low coherence (or high discontinuity) voxel points strictly, 
which could result in an extremely jagged surface. Normally, better results are 

25 obtained when the contribution from the internal force to the weighted sum in Eq.(l) 
is approximately equal to the contribution from the external force. However, the 
normalizations of the two forces may be different. Accordingly X is used as a scaling 
factor to achieve this balance. For example, when the two terms on the right-hand 
side of Eq.(3) have values in the range 0 to 255 throughout the seismic data volume, a 
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value of X in the range 0.9 to 0.95 has been found to provide a reasonable compromise 
between the internal and external forces in Eq.(l). 

Referring to step 9 of Fig. 1, values 2/ v +l and 2/ A +l represent the vertical and 
horizontal lengths respectively of a window from which internal force (curvature) and 
5 external force (discontinuity value) are computed at a fault surface point s. A 
computed energy at a surface point s will be very sensitive to noises in discontinuity 
values if the window size is too small such as (/ V =/ A =l). A good reference window 
size is the size of a window where a significant curvature can be recognized visually. 
In general, / v =5 and l k =5 (window size of 1 1 by 1 1 voxel units) provide good results. 
10 The values of l v and l h should preferably be between 1 and the maximum resolution 
the operator can reasonably expect to obtain. 

The parameter d is the maximum search distance from a point s in one search 
iteration. For a large value of d, a deformable surface is allowed to travel far distances 
in search of a desired fault surface. Too large a value of d (such as, d > 30) invites the 
15 risk of creating a surface that is quite different, in shapes and locations, from the 
initial surface. A small value of d will limit the search region. A value in the preferred 
range of 7 and 20 generally provides good results. 

Values 2r v +l and 2r A +l represent the vertical and horizontal lengths of a 
window defining the neighbors around a point s = (x,y,z). The movement of a point s, 
20 d" , is computed as a weighted average of all the values of d] in this window. This 
smoothing method prevents neighboring surface points from moving in opposite 
directions. Typical values of r v =5 and r h =5 (window size of 11 by 11 voxel units) 
provide good results. The values of r v and r h should preferably be between 1 and the 
maximum resolution the operator can reasonably expect to obtain. 

25 The step size a (0 <a <1) is introduced to stabilize the iterative minimization 

process. A small value of a , such as 0.1, guarantees high quality minimization at the 
cost of high computation time. An extremely large value of a , such as 1, makes the 
process converge fast to the minimum energy state at the risk of reaching an unstable 
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(oscillating minimum energy state) or an incorrect minimum energy state. A value of 
0.5 typically provides good results. 

This rapid fault interpretation method can be adjusted to work on a wide range 
of data quality. In good data, where faults are nicely imaged within the 
5 coherency/discontinuity volumes, the algorithm discussed above is used preferably 
with stiffness values in the range 0.8 < X< 1 and the best fit to the coherency 
/discontinuity data is obtained after five to ten iterations. The quality of the solution is 
a function of how well the input fault sticks describe the overall geometry of the fault 
surface and not how densely placed they are along the fault. 

10 Algorithms can also be used in bad data where coherency/discontinuity are 

dominated by noise or by contributions from stratigraphy rather than faults. In such a 
case, the algorithm discussed above uses a stiffness value of X = 1 to generate a dense 
fault interpretation independently of the discontinuity data but honoring closely the 
interpreted input sticks. 

15 Examples 

A performance example of 20 percent and 30 percent noise is provided below. 
A synthetic discontinuity data volume of 100x100x100 was created with a target fault 
surface y' and an initial fault surface jv' described by the equations: 

y = 70 - 0.2x - 0.2z + 5 ( sin(2nx/100) ) (4a) 
20 y=75-0.2x-0.2z (4b) 

Fig. 5 illustrates a curved target fault surface 51 and a flat initial fault surface 
53 without noise representing Eq. 4a and Eq. 4b respectively. Fig. 6A illustrates a 
horizontal slice of the discontinuity volume where the target fault surface 61 is 
25 recognized to contain high discontinuity and the initial fault surface 63 is shown as a 
solid black line. The target fault surface 61 has a mean discontinuity value of 30 and 
the whole discontinuity volume was corrupted by noises that are uniformly distributed 
between ± 5.1 (20% noise). Default parameter values, X =0.9, / v =5, / A =5, tf=7, 
a =o.5, r y =5, and r A =5, were used to move the initial fault surface 63 by minimizing 
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the energy function. In this example, 40 iterations were required to latch the initial 
surface to the final fault surface 65 as shown in Fig. 6B. 

To further test its robustness to noise, the discontinuity voxel volume was 
corrupted by noise that was uniformly distributed between ± 7.65 (30% noise). Fig. 
5 7A illustrates target surface 71 and the initial fault surface 73, where the target fault 
surface 71 is almost unrecognizable due to severe noises. In this case, a little stiffer 
surface (A =0.92) and larger windows for curvature computation and smoothing 
(/ v =5, / A =7, r v =7, and r h =5) were used to deal with severe noises in discontinuity 
values. Fig. 7B, Fig 7C and Fig. 7D illustrate horizontal slices of constructed 
10 deformable fault surfaces 75 at different stages of energy minimization (iterations of 
10, 30 and 49 respectively). 

Figs. 8A through 8D illustrate a case of 20% discontinuity noise and a 50% 
missing fault section and an initial fault surface 83. In this study, 50% of the target 
fault y* 81 in Eq. (4a) was erased in a random manner and the discontinuity data 

15 volume was corrupted by noises that is uniformly distributed between ± 5.1 (20% 
noise). As shown in Fig. 8A, the target fault surface 81 is discontinuous, but has some 
visible structural features. The same parameter set as in the 30% noise case was used. 
Fig. 8B, Fig. 8C and Fig. 8D shows a horizontal slice illustrating a constructed 
deformable fault surface 85 from three-dimensional voxel volume data at different 

20 stages of the energy minimization (iterations 5, 10 and 47 respectively). 

Examples from three-dimensional seismic discontinuity data 

Three-dimensional seismic discontinuity data from the Gulf of Mexico was 
utilized to illustrate an application of the method to actual data. Fig. 9A shows a 
25 sparse set of interpreted vertical sticks 91 displayed in three-dimensions with a time 
slice of discontinuity. The fault signature (dark line representing high discontinuity) 
93 appears clearly as a curved line on the discontinuity time slice. In Fig. 9B, the 
algorithm (or software) will generate an initial fault surface 95 by interpolating 
between the fault sticks. At this stage the initial fault surface does not fit the 
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discontinuity data. In Fig. 9C, the initial fault surface 95 is reconstructed to fit the 
discontinuity data. A stiffness value X =0.8 was used and the energy minimization 
was completed after 15 iterations. 

As shown in Fig. 9D, the algorithm outputted the interpreted fault surface 95 
5 as a set of fault sticks 97 (Fig. 9D) that can be edited if required by the seismic 
interpreter. The whole process from start to finish took about 10 seconds for the fault 
shown in Figs. 9A, 9B, 9C and 9D. The process is interactive and can be repeated to 
correct imperfections along the fault. Families of interpreted faults can also be 
processed in bulk. 
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CLAIMS 

What is claimed is: 

1. A method to create a fault surface from a three-dimensional seismic data 
volume, comprising: 

5 constructing an initial fault surface in three-dimensions, said surface 

containing at least two fault sticks, the fault sticks being from the same fault but from 
different slices of the seismic data volume, each fault stick being defined by at least 
two fault nodes; and 

reconstructing the initial fault surface to fit the three-dimensional seismic data, 
1 0 said reconstruction using an iterative evolution of a deformable surface model of the 
fault surface, said evolution being based on smoothness of the fault surface and a 
fault-indicating parameter of each location on the fault surface. 

2. The method of claim 1, wherein the fault nodes are obtained by manual 
interpretation of the data volume. 

15 3. The method of claim 1, wherein the data volume slices comprise a vertical 
slice and an horizontal slice. 

4. The method of claim 1, wherein the seismic data comprise seismic amplitude 
discontinuity data, and the fault-indicating parameter is the seismic amplitude 
discontinuity. 

20 5. The method of claim 4, wherein the iteration of the deformable surface model 
uses local minimization of an energy function of such surface, said function having as 
variables the curvature of the fault surface and the discontinuity in seismic amplitudes, 
each variable being itself a function of spatial location on the fault surface, and 
wherein the iteration is terminated when the change in the surface energy relative to 

25 the preceding iteration is less than a pre-determined value. 

6. The method of claim 5, wherein the energy function is a weighted sum of an 
internal force and an external force, said internal force being a function of curvature of 
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the fault surface and said external force being a function of seismic amplitude 
discontinuity on the fault surface. 

7. The method of claim 1, wherein the seismic data comprise seismic amplitude 
coherency data, and the fault-indicating parameter is the seismic amplitude coherency. 

5 8. The method of claim 7, wherein the iteration of the deformable surface model 
uses local minimization of an energy function of such surface, said function having as 
variables the curvature of the fault surface and the coherency in seismic amplitudes, 
each variable being itself a function of spatial location on the fault surface, and 
wherein the iteration is terminated when the change in the surface energy relative to 
10 the preceding iteration is less than a pre-determined value. 

9. The method of claim 8, wherein the energy function is a weighted sum of an 
internal force and an external force, said internal force being a function of curvature of 
the fault surface and said external force being a function of seismic amplitude 
coherency on the fault surface. 

15 10. The method of claim 6 or claim 9, wherein the weighted sum uses a weighting 
factor X for the internal force and a weighting factor (1 - X) for the external force, X 
being selected to provide substantially equal contributions to the energy function from 
the internal and external forces, except for noisy data where X is chosen substantially 
equal to 1 . 

20 11. The method of claim 1, wherein a fault surface is assumed to be a single- 
valued function of two orthogonal geographic coordinates, and all voxels in the fault 
surface are allowed to move only in the third orthogonal direction in the course of 
iterative deforming of the surface model, said single- valued function being in the form 
of x = f(y,z) if the angle between a fault surface and the y-z plane is less than 45 

25 degrees, and y = f (x,z) if the angle between a fault surface and the y-z plane is greater 
than 45 degrees, the coordinate z being in the vertical direction. 

12. The method of claim 1, wherein constructing an initial fault surface comprises 
the steps of: 
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(a) generating a pre-selected density of fault points along each fault stick by 
performing linear interpolation between fault nodes; 

(b) identifying the extent of the boundaries of the fault surface by using the 
interpolated fault points from step (a); and 

5 (c) constructing an initial fault surface within said boundaries by using 

bilinear interpolation of the interpolated fault points from step (a). 

13. The method of claim 1, wherein the reconstructed fault surface is defined by a 
discrete set of fault sticks. 

14. The method of claim 1, wherein the seismic data comprise seismic amplitude 
10 discontinuity and coherency data, and the fault-indicating parameter is a combination 

of the seismic amplitude discontinuity and the seismic amplitude coherency. 



G 



G 



WO 2004/038654 



1/5 



Seismic discontinuity 
data volume 



PCT/US2003/029424 

1U/ 527 4 15 



Seismic amplitude 
data volume 



Manually picked fault 
surface points 



Initial fault surface generation 
through interpolation 



X 



9 



Reconstruction of a fault 
surface to fit the seismic 
discontinuity data 



FIG. 1 




r23 



FIG. 2B 



FIG. 2 A 



( 



b 11 MAR 2005 



■^vMfcmh 3LAW8( (usptc, 




FIG. 4 A 



FIG. 4B 



( 

tUM tm a rw iir- 1 k) J J MAR 2005 



c . c 

%0: 2004/038654 PCT/US2003/029424 

AO/ 527 4 15 

3/5 





FIG. 6 A 



FIG. 6B 



( 

STDS Rafil FCT/PTQ 



11 MAR 2005 



IK (USPTC) 



c 

W&2004/038654 



C 10/527415 

PCT/US2003/029424 





( 

mm Kec'o mm® i \ mar 2005 



THIS -rmt nuhMK 



c 

SHTORec'dPE., • j j 




This Page is Inserted by IFW Indexing and Scanning 
Operations and is not part of the Official Record 

BEST AVAILABLE IMAGES 

Defective images within this document are accurate representations of the original 
documents submitted by the applicant. 

Defects in the images include but are not limited to the items checked: 

□ BLACK BORDERS 

□ IMAGE CUT OFF AT TOP, BOTTOM OR SIDES 

□ FADED TEXT OR DRAWING 

Of BLURRED OR ILLEGIBLE TEXT OR DRAWING 
0 SKEWED/SLANTED IMAGES 

□ COLOR OR BLACK AND WHITE PHOTOGRAPHS 

□ GRAY SCALE DOCUMENTS 

□ LINES OR MARKS ON ORIGINAL DOCUMENT 
□TrEFERENCE(S) OR EXHIBIT(S) SUBMITTED ARE POOR QUALITY 

□ OTHER: 

IMAGES ARE BEST AVAILABLE COPY. 
As rescanning these documents will not correct the image 
problems checked, please do not report these problems to 
the IFW Image Problem Mailbox. 



■ riiS-PAGE BLANK (usptg, 



